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Abstract 

The flow field is simulated on the surface of a Butler wing in a uniform 

l 

stream. Results are presented for Mach number 3.5 and Reynolds number of 
2,000,000. The simulation is done by integrating the viscous Navier-Stokes 
equations. These equations govern the unsteady, viscous, compressible and heat 
conducting flow of an ideal gas. The equations are written in curvilinear 
coordinates so that the wing surface is represented accurately. 0-type and H- 
type grids have been used for this study, and results are compared. The 
governing equations are solved by the MacCormack time-split method, and the 
I results are compared with other theoretical and experimental results. The 

j codes are written in FORTRAN, vectorized and currently run on the CDC Vector 

I Processing System (VPS-32) computer. 
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Nomenclature 


c v 

e 

E 

F, G, H 
J 
k 
K 
L 
M 
n 
P 

Pr 

V ^y* ^2 

R 

R e 

T 

T. 

T w 

T r 

t 

s o 

0 

u,v,w 


x.y.z 


a 


specific heat at constant volume 
internal energy per unit volume 
total energy per unit volume 
vector fluxes for coordinates directions 
Jacobin matrix 

coefficient of bulk viscosity 
coefficient of thermal conductivity 
length of the wing 
Mach number 

magnitude of normal vector on any surface 
static pressure 
prandtl number 

components of heat conduct vectors 
universal gas constant 
Reynolds number 
static temperature 
free stream temperature 
wall temperature 

reference temperature for Sutherland viscosity law 
time 

constant for Sutherland viscosity law 
vector of state variables 

velocity components in the physical coordinates 
contravariant velocity in the computational coordinates 
coordinates for the physical domain 
coefficient for pressure damping 


v 


Y 

5.ti,C 

n 

*V 

P* 

P 

T 

T , X , T 

xx’ xy’ xz 
T yy’ T yz* T zz 

WV 6 ij 


ratio of specific heat 

coordinates for the computational domain 

molecular viscosity 

reference viscosity in Sutherland viscosity law 
bulk viscosity 
density 
stress tensor 

elements of the stress tensor 

components of viscous dissipation functions 

Kronecker delta function 


subscript 


w 


solid walls 
free stream value 


superscript 

n time level 


indices 

i,j,k point indices 


operators 
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o 

L5.Ln.LC 

a 


gradient 
inner product 

finite difference operators 
partial differentiations 
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Introduction 


The Butler wing is a delta wing which was proposed by D. S. Butler 1. The 
planform of the body is an isosceles triangle, and the leading edges of the 
wing lay along the Mach lines of the unperturbed stream. The first 20% of the 
wing is conical and the last 80% of the wing has elliptical cross sections 
with increasing eccentricity along the x-axis. At the trailing edge, the 
elliptical cross section has infinite eccentricity and this last cross section 
is a straight line. D. S. Butler has compared the experimental results for 
surface pressure with the theoretical results using the slender-body theory 
approximation to simplify the inviscid equations of motion*. Walkden and Caine 
estimated the pressure on the surface of a Butler wing at zero incident in a 
steady uniform stream by numerically integrating the two semi-characteristic 
forms of the equations governing inviscid supersonic flow of an ideal gas with 
constant specific heat^. Squire obtained experimental results for a Butler 
wing varying Mach number and angle of attack^. Inviscid equations have been 
used in all previous analytical and numerical investigations. 

In order to study the flow around a Butler wing, the Navier-Stokes 
equations are numerically solved. The equations are unsteady, compressible, 
viscous and three-dimensional. The time dependency of the governing equations 
allows the solution to progress naturally from an arbitrary initial guess to 
an asymptotic steady state, if one exists. The equations are transformed from 
physical coordinates to computational coordinates which allows solutions to be 
computed in a rectangular domain. The equations are solved by the MacCormack 
time-split technique'** ^ which is vectorized and programmed to run on the CDC 
VPS-32 (CYBER 205) computer. The codes are written in 32-bit (half-word) 
FORTRAN. 


The Butler wing is symmetric about (x-z) and (x-y) planes. This permits 
the use of one quarter of the entire physical domain for flow field 
computation with zero angle of attack (Fig. 1). However, if the angle of 
attack is greater than zero, then half of the physical domain should be 
considered. 

Grid generation is the first step which should be considered in obtaining 
flow field solutions over any configuration. Due to the data base management 
of the present program, it is necessary to map entire physical domain into a 
rectangular parallelepiped. Among the grid types, selection of an 0-type grid 
for cross sections in the x-coordinate direction would produce a point 
singularity at the nose tip and a line singularity along the trailing edges. 
Nevertheless, an 0-type grid maps the solid boundary onto an entire face of 
the parallelepiped. It is also possible to generate a highly orthogonal grid 
in the regions where there are relative high curvatures. However, a H-type 
grid does not map the solid boundary onto an entire face of the computational 
box. This creates a potential problem in updating the boundary conditions near 
the leading edges of the wing and also the grid in some regions could be 
highly skewed. But, there are no singularities in the grid. Both types of 
grid have been used in this study, and results are compared with other 
numerical, analytical and experimental results. 
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Governing Equations 


The governing equations for a thermal fluid system are the conservation 
of mass, momentum and energy. These equations are developed for an arbitrary 
region under the assumption that the system is a continuum. Equations of 

motion for viscous, compressible, unsteady and heat conducting flow can be 
written as: 


Continuity: 

|| + v*(pu) » 0, 

(la) 

Momentum: 

^ ♦ 7.<p55 - ?) - 0, 

(lb) 

Energy: 

+ V*(Eu + q - u * t) = 0, 

(lc) 


where E is the total energy per unit volume given by E = p (e + u»u/2) and e 
is the internal energy per unit volume. Equations lb and lc can be simplified 
by assuming that the stress at a point is linearly dependent on the rate of 
strain (deformation) of the fluid (Newtonian fluid). 


du.j du 1 . du 

T i j = * p6 ij + p( 3xT + *x? } + 6 i j 


(2a) 


where 6.. is Kronecker delta function, and p' is the second coefficient of 

• J 

viscosity which is related to the coefficient of bulk viscosity (x) by the 
expression k= 2 p/3+ p'. The contribution of ic can be neglected if the 
pressure in a fluid is not changed abruptly during its expansion or 
contraction. Under this assumption, the stress tensor can be related to the 
pressure and velocity components as: 



For an isotropic system, the heat flux in Eq. (lc) can be expressed in terms 
of temperature gradient (Fourier law of heat conduction) as: 

q = - KVT (3) 


where K is the coefficient of thermal conductivity. A common approximation 
used for viscosity is based on the kinetic theory of gases using an idealized 
intermolecular-forces potential; the relation is 


where 



3/2 T r + S o 

T+T' 


S rt = 198 *R; 
o 

u * 0.1716 np 
r 


(4) 


The coefficients of thermal conductivity K can be determined from Prandtl 
number as 


where C y is the specific heat at constant volume and y is the ratio of 
specific heats. 

It is essential to have a supplementary relation to close the system of 
equations (la-lc). By neglecting the inter-molecular forces (thermally perfect 
system), thermodynamic properties can be described as: 


P = PRT, 


( 6 ) 


where R is the universal gas constant. Thermally perfect gas assumption 
permits to express the internal energy as a function of temperature only i.e., 
e=e(T). In addition, assumption of calorically perfect gas Ce (0) =03 allows the 
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following relation: 


e = C y T. (7) 

A substitution of Eq. (6) into Eq. (7) results in 

p = pe (y-1). (8) 

The equations of motion are in conservative form. For simplicity, these 
equations can be written in a compact vector form as 


where 




■ 


p 
pu 
pv 
pw 
E 

V j 


> 


dU . dF . dG . 5H n 

3I + 3x + 3y + 3z =0 ’ 


r 


pu 


puu - x xx + P 


. F - \ 


puv - T, 


xy 


puw - T, 


X2 




Eu + q - ♦ + pu 

V x X y 


(9) 


G = 


/ 





pv 



pw 


puv 

X 

1 


puw - 

T zx 

pw 

- T + P , 

yy K 

r { 

pvw - 

T zy 

pvw 

N 

1 

V 

pww - 

T + D 

ZZ V 


+ \ - * y + 

Ew + 
V 

q z - 4> z + pw 

J 




For the sake of generality, the governing equations are transformed from a 
physical domain into a computational domain as 


5 


r 



r dF dG dH > 

W Sp w = 


o. 


( 10 ) 


The transformation coefficients can be computed from a functional 
relation between the computational coordinates and the physical coordinates. 

5 * S(x,y,z), t) = n(x,y,z), C = C(x,y,z). (11) 

If Eq. (11) were known, the transformation coefficient could have be computed 
by direct differentiation. If not, after some algebraic manipulations, the 
transformation coefficients can computed by 

«x= (y „ 2 c ' y c S’ ' I 0 ' 1 !- 

5 y = '1% " x ; 2 t,* I I J I" 

x C y n ) ' I 0 * 1 !' 

, - 1 , 1121 

h x (y^ z^-y^z^) / |J | , 

■Hy= (x^Z^-X^Z^) / |j |, 

\ = “ (x c y C “ x c V / 

Cx = (y5 zt| - yn zO / |J I* 
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S ' " (X C Z >1 " X T| * l J I* 


? 2 ’ (x 5 V ? > ' I 0 ' 1 !’ 


where [d _1 ] is defined as 


J’ 1 ! * x 5 (y „ \-H V ‘ % iy l z c' y C z 5 ) 


* X C (y 5 %' ’V? 1 ' 


(13) 


In the present case, the grid planes are perpendicular to the x-coordinate. 
Consequently, physical coordinates can be written as 
x = xU) , 

y = yU.n.c), (14) 

z = z(5,T|,C). 


This reduces the transformation coefficients from nine to five none zero 
elements therefore reducing the memory requirements. 
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Method of Solution 


A time marching method is used to compute the solution so that the 
possible transient features can be readily captured. This method is an 
explicit second-order accurate time-split predictor-corrector algorithm The 
governing equation, Eq. (10) is discretized in computational directions. In a 
compact form, it can be expressed as 


where 


- n+1 
U 

i jk 


= CL T) (At T| )]CL c (At c )]CL 5 (At 6 )]CL c (At c )]CL T) (At T) )]. (15) 


A S = At c - 7 A V 


and L , L and L , are the operators in £, n and C directions, respectively. 
A time step is completed in this algorithm with the application of each 
operator applied symmetrically about the middle operator. For example, 
operator can be defined as 

out 

L-(AtJ = 0 , (16) 

5 5 1 jk 


for predictor step: 


- in 

0 = u 

i jk i jk 


At C 

TC 


[ (F . 



a? 

a* 


* (5 i - 5 i-i>f + (n i - W # ] i 




jk 


(17) 


for corrector step: 


out 

U 

i jk 



i n At- 



ae 

*5x 


(S i+l ‘ V |y + (R i+l “ R i } U- 1 . (18) 

i jk 
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The solution is stable if the time step of each operator does not exceed the 
allowable step size for that operator. The finite difference scheme is 
consistent if the sum of the time steps for each operators are equal. The 
solution is second order accurate if operators are applied symmetrically. 

This method has a time-step stability limit, but there is no rigorous 
stability analysis available. A commonly used conservative time-step is 


At < min 


Ax Ay Az 





( 19 ) 


where c is the local speed of sound. 

In the supersonic region, there exits a large gradient which requires a 
very fine mesh to resolve it. Most center difference methods admit a solution 
which has sawtooth or plus-minus waves with the shortest wave length that the 
mesh can support. In the case of a nonlinear problem, these short waves 
interact, vanish, and reappear again as distorted long waves or 
oscillations. These oscillations eventually blow up the solution, if they are 
not resolved. The oscillations of "low frequency" can be suppressed by adding 
a fourth order damping. A common damping used is the pressure dampening. This 
can be expressed in computational coordinates as 


- a 


A 


At. 


3d V A + c d 2 p dU 


% = 1 , 2,3 


( 20 ) 


where VA is the contravariant velocity. 


Initial and Boundary Conditions 


In computational fluid dynamics the initial conditions usually correspond 
to a real initial situation for a transient problem, or a rough guess for a 
steady state problem. In practice, initial conditions are obtained from 
experiments, empirical relations, approximate theories or previous 
computational results. An improper initial guess may result in generating 

unrealistic strong transient waves which propagate through the computational 
region dominating the flow field and eventually lead to a solution failure. An 

important requirement for the initial conditions is that they should be 

physically as close as possible to the actual nature of the flow field in the 
region under study. This will minimize the number of iterations required for 
convergence. An attractive approach is to initialize the entire flow field 
with a crude and simple guess (e.g., free stream condition). During the course 
of the computation, both body and upstream boundary conditions are changed in 
a gradual manner to their final values over a prescribed number of iterations. 
In the present study, this technique is applied in only one step which is 
equivalent to impulsive initial conditions. 

It Is equally important to implement a realistic, accurate and stable 
method to determine boundary conditions. The application of certain conditions 
may cause numerical instability even though the flow is physically stable. 
Most of the boundary conditions currently implemented are drawn mainly upon 
intuitions, simple analytical expressions, wind tunnel experiments and 
computational experimentations. In the selection of boundary conditions, 

consideration should be given to the following criteria: convergence, 

stability, computer time and above all the physical justification. 

For the present case, there are five different boundary conditions. They 
are upstream, downstream, lateral, top, and solid boundary. For the case of H- 


grid, the upstream boundary is located at six grid point spacings ahead of the 
nose of the wing. Following undisturbed free stream conditions are assumed 
for this boundary: 


U 


upstream 


= 0 


oo 


( 21 ) 


For the case of O-grid, the upstream boundary is set at 5 % of the body length 
to avoid the singularities. The conical assumption has been made for this 
boundary®. Flow is said to be conical if the physical conditions such as 
pressure and velocity do not vary with position along any ray through a point, 
which is referred to as the vertex. For this case, the viscous-conical 
solutions are obtained for a cone at proper angle of attack. This is done by 
creating a conical grid, which is straight lines (rays) from the vertex and 
integrating the conical Navier-Stokes equations. This solution is valid 
provided the body is sufficiently slender. 

A zero gradient in the y-direction (parallel to the primary direction of 
flow) is assumed for the downstream boundary, i.e.. 


au 

* 5 * 


= 0 . 


( 22 ) 


The lateral boundaries are located far enough away form the body to avoid any 
influence on the interaction region. Presently, a zero gradient in the z- 
direction is assumed for these boundaries, i.e.. 


aO 

■§z 


Top 


= 0 . 


( 23 ) 


The walls are assumed to be impermeable and no-slip boundary conditions are 
applied, therefore, all velocity components are assumed to be zero. Similarly, 
the wall is assumed to have a constant temperature T w . A zero normal pressure 
gradient is assumed for the solid surface, i.e.. 
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dP 

3n 


Solid 


= 0 . 


(24) 


This appears to be a boundary-layer approximation (i.e., a zero normal 
pressure gradient). It is, however, a much milder approximation, since 
constant pressure is not applied through the boundary layer but over one grid 
line in the boundary layer. This approximation has yielded stable computation 
for both the non-separated and separated boundary layers^. For general 
curviliear coordinates, Eq. (24) can be expressed as 


\ + **y ’V + | + 

(n 2 + n 2 + n 2 ) P + 

x y z n 

( ?x \ + c y "y + C z *) P ? ' °‘ 


(25) 


1 This equation is satisfactory even for skew grids near solid boundaries. 

I 
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Applications to a Butler Wing 


As mentioned in the introduction, the Butler wing is a delta wing which 
was proposed by 0. S. Butler 1 . The Butler wing is analogous to an elliptic 

j piston which the ratio of major and minor axes changes while disturbing the 

fluid. The plan form of the body (Fig. la) is an isosceles triangle, and 
leading edges of the wing are laid along the Mach lines of the undisturbed 
stream. For the first 20% of the length, the body is a right circular cone. 
The reminder of the body has elliptical sections which become more eccentric 
as the sharp trailing edges are approached (Fig. lb). The semi major and minor 
axes are given by: 

Major axis (semi -span) = 

Minor axis (thickness on = 
centerline) 

I 

I 

! where £ - M^-l 

The model is 0.8 ft. (0.2438 m) long, and the geometry has been generated 

i 

i for Mach number of 3.5. That is the semi-apex angle of the plan form and the 

initial conical nose is SIN -1 (l./3.5) =16.602°. 

i This wing is symmetric about (xz) and (xy) planes. For zero angle of 

attack, this permits us to use only one quarter of the entire physical domain 
(Figs, la-lc) which is extremely advantageous from computational viewpoint. 
However, if the angle of attack is non-zero then half of the physical domain 
should be considered. 

Some specific flow field results have been obtained for the Butler wing 
and these are discussed in the next section. 


x 

? 

x 

? 


X [1 


? 


- ( 


x- 0.2L \ -\ 

of' j 


(26a) 

0<=x<=l 

(26b) 

0<=x<=0.2L 


(26c) 

0.2L<=x<=L 
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Results and Discussion 


In the previous experiments'*, the surface oil-flow patterns over the wing 
at various angles of attack and at Mach number 3.5 showed no signs of 

transition and the nature of the oil streaklines was typical of a laminar 

flow. Therefore, results are obtained only for laminar flow over a Butler wing 
at Mach number of 3.5, Reynolds number of 2xl06/ft (6.56xl06/m), free static 
temperature of 390°R (216.67®K), wall temperature of 1092®R (606.67°K), length 
of 0.80 ft (0.2438 m), and at zero and ten degree angles of attack. In this 
study, a two boundary grid generation (TBGG) technique 4 is used, this method 
is essentially an algebraic method. The application of the TBGG method 
requires that the entire body be sliced into different cross-sections. These 

cross-sections are obtained in the stream-wise direction by analytical 

descriptions of the wing surface (Eqs. 26a-26c). Then, two types of grid 
generated for this wing the H-type and 0-type. Results of both cases are 
compared and discussed. 

H-Type Grid 

In this case, the entire flow field is sliced into fifty-five stations in 
the stream-wise direction, and each station has 64x36 grid points (Fig. 2). 
There is a total of 126,720 grid points which take 2.8 million 32-bit words of 
primary memory (16 variables). The required computational time is 1.9x10-5 
sec/grid point/i teration (2.5 sec/iteration) which is typical for the CYBER 
205 with two pipes. Results are obtained for zero angle of attack. The 
computed pressures are plotted in Figs. 3a-3c. The pressure coefficient along 
the center line is shown in Fig. 3a. The results are compared with available 
experimental and numerical results (Refs. 1,3,8). The results on the center 
line are in excellent agreement with the experimental and previously obtained 
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numerical results (Refs. 1,3,8). At 41.66% and 68.33% chordwise position, the 
pressure ratios are plotted against the conical spanwise coordinates y/xtan 
(<t>). Figs. (3b and 3c). They are in good agreement with experimental and 
numerical results. However, There are some discrepancies in the results 
between 30° and 60®. This is probably due to the fact that grid lines are 
notorthogonal near those regions and is a direct consequence of the H-type 
grid. 

0-Type Grid 

For this case, the physical domain is limited to 5% to 95% of the wing. 
This is done to avoid any singularities. The conical Navier-Stokes solutions 
are forced for the upstream boundary which is located at 5% of the wing. This 
solution is obtained by the integration of the Navier-Stokes equations 6 for a 
conical grid with proper boundary conditions. The wing is sliced into forty- 
one stations in the stream-wise direction, and each station has 41x127 grid 
points (Fig 4). There is a total of 213,487 grid points which take 4.7 million 
32-bit words of primary memory. Results are obtained for zero and ten degree 
angle of attack. 

Results for zero angle of attack are compared with results obtained with 
the H-type grid, from the experiments and other numerical results (Refs. 
1,3,8). The computed pressures are plotted in Figs. 5a-5c. The pressure 
coefficient along the center line (Fig. 5a) is in good agreement with other 
numerical and experimental results. Nevertheless, there is some discrepancy 
near the nose region. This may be due to the fact that the upstream solutions 
are based on the conical solutions. But, solutions match exactly with results 
from H-type grid, this is because grid topology near the center line is the 
same for both grid types. At 41.67% and 68.33% chordwise position, the 
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pressure ratios are plotted against the conical spanwise coordinates y/x tan 
(<}>). They are in excellent agreement with experimental and numerical results 
(Fig. 5b and 5c). In addition they are much closer to the experimental results 
compared with results from H-type grid. This is due to good grid orthogonality 
in the case of 0-type grid. On the thick sections near the nose the pressure 
is highest on the centerline and falls toward the leading edge. Figure 5d 
shows the cross-flow velocity at 5%, 23%, 412, 59%, 77 % and 95%. 

Results for ten degrees angle of attack are compared with experimental 
results. The computed pressures are plotted in Figs. 6a-6d. At 17%, 30%, 50% 
and 70% chordwise position, the pressure ratios are plotted against the 

conical spanwise coordinates y/x tan ( 4> ) . They are in good agreement with 
experimental and numerical results (Fig. 6a-6d). On the thick sections near 
the nose the pressure is highest on the centerline and falls toward the 

leading edge whereas near the trailing edge the spanwise distribution is more 
'wing like' with the maximum pressure at the leading edge. The changeover is 
shown by the pressure peaks in the pressure distributions at x/c=0.5 and 0.7 
at 10° angle of attack. There are some discrepancy near x/c=30-50; this may be 
due to the fact that Squire 3 has not used the exact model of the Butler wing. 
In order to mount the model in the wind tunnel, the lower surface was 

distorted to Include a sting support. Figure 6e shows the cross-velocity at 
5%, 23%, 41%, 59%, 77% and 95%. These figures show a weak cross-wise 
separations at suction side which is confined to the body. At 59%, the cross 
flow has separated but a well defined vortex is not visible. Squire® has 

performed a series of tests to investigate the effects of thickness on the 
longitudinal characteristic of a delta wing of aspect ratio. The tests on the 
thick symmetrical delta wings have confirmed that the lift curve slope 
decreases as thickness is increased. This loss of lift is associated with a 
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weaker vortex system giving less nonlinear lift. Squire 3 also observed a pair 
vortices at the trailing edges, but there was no sign of any spanwise flow 
outboard of these vortices. 
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Concluding Remarks 


General formulations are presented to investigate the flow field over 
complex configurations for high speed free stream conditions. An advanced 
algebraic method is used to generate grids around these configurations. The 
computational procedure developed is applied to investigate the flow field 
over a Butler wing. Illustrative results obtained for specified free stream 
conditions compare very well with available experimental and numerical 
results. Results are obtained only for laminar flow over a Butler wing at Mach 
number of 3.5, Reynolds number of 2xl06/ft (6.56xl06/m), free static 
temperature of 390°R (216.67°K), wall temperature of 1092 # R (606.67°K), length 
of 0.80 ft (0.2438 m) and the wing is at zero and ten degree angles of attack. 
Two types of grid have been generated for this wing; H-type and 0-type. 
Results of both cases are compared and discussed. 
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Figure 6a. Pressure Coeffient 
(X/C 17% , Ten Degree Angle of Attaek) 
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Figure 6b. Preseure Coeffient 
(X/C 30% , Ten Degree Angle of Attaek) 





Figure 6e. Preseure Coeffient 


Figure 6d. Pressure Coeffient 


(X/C 60% , Ten Degree Angle of Attack) 


(X/C 70% , Ten Degree Angle of Attack) 
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x/c=59% 
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Figure 6e. Pressure Coeffient (X/C Ten Degree Angle of Attack) 


